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Abstract. Unitary transformations can allow one to study open quantum systems in 
situations for which standard, weak-coupling type approximations are not valid. Here, 
we develop an extension of the variational (polaron) transformation approach to open 
system dynamics, to apply to arbitrarily large exciton transport networks with local 
environments. After deriving a time-local master equation in the transformed frame, 
we go on to compare the population dynamics predicted using our technique with other 
established master equations. The variational frame dynamics are found to agree with 
both weak coupling and full polaron master equations in their respective regions of 
validity. In parameter regimes considered difficult for these methods, the dynamics 
predicted by our technique are found to interpolate between the two. The variational 
method thus gives insight, across a broad range of parameters, into the competition 
between coherent and incoherent processes in determining the dynamical behaviour of 
energy transfer networks. 



PACS numbers: 31.15.xt, 31.15.xp, 03.65.Yz 
1. Introduction 

The theory of open systems is necessary to describe any quantum system in contact with 
an uncontrollable and non-negligible environment. In problems of energy transport one 
is often interested in a regime where the environment, which consists of a very large 
number of degrees of freedom, is highly influential. We are here concerned with the 
dynamics of electronic excitations across some discrete network of molecules, in which 
the environment can play a key role - and this is exemplified by recent observations in 
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various microscopic biological systems. Examples include: Molecular complexes involved 
in photosynthesis, such as the Fenna- Matthews- Olsen (FMO) complex in green sulphur 
bacteria (HEHSHUElEllTllHlini HDl HI] and the light harvesting photosystems in green 
plants [12]; the electron transport chain in Respiratory Complex I ^13j; the donor- 
bridge-acceptor model of olfaction [H] and certain models of magnetoreception in birds 
[T5| IT6] . Similar energy transport models have also been used extensively outside of 
biology - applications range from the dynamics of coupled quantum dots [T71 HH [El [20] 
in sohd state physics, to those of impurities in lattice Bose-Einstein condensates (BECs) 

Several techniques have been developed over the years to calculate the explicit time 
domain dynamics of open quantum systems. Some are numerically exact, meaning that 
given sufficient computational resources, they will converge to the correct dynamics 
under some well-controlled approximations. Such techniques include the path integral 
[22| |23| [21], hierarchical equations of motion (HEOM) [U [251 [21] and density matrix 
renormalization group (DMRG) [27] methods. Though powerful, these approaches 
typically place restrictions on the kind of system that can be modelled, and they may 
also scale badly (in terms of computing resources) with the size and complexity of said 
system. Often, as will be the case in this article, more numerically tractable, though 
approximate methods are used, such as those based on master equations [281 [211 [30] . 
This technique provides an equation of motion for the reduced density matrix of the 
system in question without having to track the full evolution of the environment, though 
normally involves some kind of perturbative expansion in a small parameter, such as 
the system-environment coupling strength. A further advantage of the master equation 
approach is that it can offer insights into the mechanisms underlying the dynamics of a 
system by relating rates and energy shifts directly to microscopic parameters. However, 
the obvious drawback of many master equations is that they rely on certain Hamiltonian 
parameters being small. If this condition is not fulfilled, then the truncation of the 
perturbative expansion often leads to (potentially unphysical) results which can diverge 
wildly from the true dynamics [31] . 

In certain parameter regimes, performing unitary transformations, such as the 
polaron transformation [32l |33l [HI [THl [35l [36l [371 [33 [39l [40] , on the combined system- 
environment Hamiltonian can result in a smaller interaction energy in the transformed 
frame. The transformed system is then amenable to being modelled using a perturbative 
master equation. For example, the polaron transformation can work well over a broad 
range of parameters when the relevant environmental timescales are short compared to 
those in the system - in fact, the polaron transformation diagonalizes the Hamiltonian 
we use below when no electronic couplings are present between the sites. It is thus 
often used when the coupling between system and environment is strong. Between 
the weak-coupling and polaron regimes, however, lies a region of parameter space for 
which neither model is appropriate. In addition, the polaron transformation runs into 
problems when applied to a system with an Ohmic or sub-Ohmic environment (one for 
which the environment spectral density scales linearly or sub-linearly, respectively, at 
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low frequencies). In this case, infrared divergences arise which prevent certain master 
equations (such as the time-local form used in this work) from correctly predicting the 
dynamics in the transformed frame. 

As an extension of the standard polaron approach, the variational (polaron) 
transformation [HI |l2l HHl HU HHl SHI HT] allows one to attempt to derive a perturbative 
series which is as valid as possible (given the restricted form of the transformation) 
in all parameter regimes. This is achieved by performing an optimized, partial 
polaron displacement on each of the environmental phonon modes relevant to their 
particular mode frequency, thereby interpolating between the weak-coupling and polaron 
representations for separate modes, as well as in the final master equation. Here we build 
on previous work on the variational transformation for two-site systems |l3lllHll6] which 
is in turn based on an idea originating with Silbey & Harris [HI |42]. The major new 
contribution of this paper is the generalization of the formalism to any number of sites, 
allowing for the simulation of large networks across a range of environmental coupling 
parameters and temperatures. For comparison to other techniques, we have included 
examples of dynamics for systems in several different regimes. 

In section [2] we describe the mathematical model for which the transformation is 
valid, including some of its limitations. We then go on to discuss the form of the 
variational polaron transformation and the accompanying optimization procedure in 
section [3l In section H] we outline the master equation formalism in the variationally 
transformed frame, and in section [5] we present some example dynamics, including that 
of the FMO system. Finally, in section El we shall conclude by briefly discussing the 
various pros and cons of the method outlined herein. 

2. The transport model 

The system (S) considered in this work is that of coupled two-level systems, known 
as sites. Between them they carry exactly one excitation - for molecular networks these 
are electronic excitations which, for charge neutral systems, are called excitons. The 
latter restriction to a single excitation allows for a great reduction in the size of the 
system Hilbert space, and is sufficient to describe the behaviour of many physically and 
biologically relevant systems. For example, it is thought to be a valid approximation 
to the in vivo dynamics of the FMO complex studied in section [5] [5J. Each of the 
sites interacts with its own, independent phonon environment (E) through a linear, 
position-based coupling. The Hamiltonian for the combined system and environment is 
given by 

H = Hs + He + Hj, Hs = Y.^n\n){n\+ V;^|n)(m|, 

He = ^ Wn,k&i,A,k, Hi =Y, l^)(^l(^n,k&l,k + ^n.k^n.k), (l) 

n,k n,k 

where fe„^k is the annihilation operator for phonon mode k on site n, and \n) is the state 
of S in which only site n is excited (see figure [1] for a cartoon visualization of the above 
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Figure 1. A schematic representation of the exciton transport network. The numbered 
ellipses represent the acceptor sites of the system, each of which is coupled to its own 
environment of phonon modes (masses on springs) and, in general, to every other site 
- via the black lines. 



Hamiltonian) . 

The environmental spectrum is usually taken to be a continuum, such that the 
couplings gn,k can be described in terms of a spectral density J„(a;). The spectral 
density is a continuous function of frequency defined by 



which takes into account the density of states, dispersion relation and interaction 
mechanism with the environment. In the continuum case, a good measure of the strength 
of the system-environment coupling at each site is the reorganization energy: 



This model ignores any spatial correlations between phonon excitations at different sites, 
meaning that the Hamiltonian in [1] is not relevant for systems with strong, long-range 
correlations, such as impurities in BECs. For the case of FMO it has been claimed, 
based on detailed molecular dynamics simulations, that spatial correlations do not play 
a significant role in the exciton dynamics jUj. Since the environments at each site 
are independent, one can take the couplings to be real for this model without loss of 
generality. In the case of a global environment, however, the phases of the couplings to 
each site can encode the environmental correlations between them. 

It is generally assumed that the environment is initially in a thermal (Gibbs) state 
at temperature T = l/^ksf^): P£;(0) = e~^^-^/tr(e~^^-^). In addition, the combined 
system-environment state is assumed to be initially separable, such that there are no 
system-environment correlations: p(0) = ps(0) ®P£;(0). 




(2) 



k 




(3) 
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3. The variational polaron transformation 

The Hamiltonian in [1] describes the transport problem in a way which is intuitive and 
transparent, in that Hs and He are the Hamiltonians for the system and environment 
in isolation and the interaction term Hj is simple in form. However, one of the features 
of quantum mechanics is that the physics of a system - even a composite one - is 
invariant under unitary transformations. As a consequence, one is not restricted to a 
single way of distinguishing two subsystems. In the model outlined above, much of the 
energy associated with the interaction between system and environment is due to the 
excitation deforming the surrounding molecular structure, and hence affecting the state 
of the phonon environment. 

Applying the polaron transformation allows us to move into a reference frame where 
this back-action from system to environment is accounted for at the Hamiltonian level. 
The system is 'dressed' by the environment, and the environmental phonon modes are 
displaced in position conditional on the state of the system. The explicit form of the 
transformation is 



H 



where G = ^ |n) (n|u;„j,(/„,k6^ 
leading to a transformed Hamiltonian: 



k"n,kj) 



(4) 



?i,k 



H = Hq + Hi, 

Hq = Hs + He, Hj = Hl + Hd, 

Hs = Yl + RrM{n\ + BnBmVnrn\n){m\, 

Hl = J2 [(9n,k - /n,k)fel,k + {9n,k " /n,k)*&n,k 

n,k 

Hd = X! ^nm\'n){m\Bnm, He = He- 



(5) 



(6) 



The interaction Hamiltonian Hi now contains two terms. One, H^, is of the same linear 
form as the interaction in the untransformed Hamiltonian, albeit with modified coupling 
strength. The other term. Ho, contains a new kind of interaction between off-diagonal 
system operators and the environmental displacement operators 



Bnm = = exp 



X! '^n,k(/n,k^|i,k fn,]J^n,k) X! '^m!k(/m,k&L,k /m,k^m,k) 



L k 

-B B 



Here, the traces {BnBm) of the displacement operators in the Bnm have been taken into 
the system Hamiltonian, and are thus treated as renormalized couplings between the 
sites. The S„'s are given by 



i3„ = tr < exp 



X!'^n,k(/",k&l,k ~ /n,k^n,k) 



L k 



PE 



exp 



1 ^ LM! eoth(/3a;„,k/2) 



n,k 
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for a thermal equilibrium environmental state pe- The site energies after transformation 
are also shifted in comparison to the original frame by a factor defined as 

k 

Usually when the polaron transformation is discussed in the literature, what is 
meant is the fully- displaced version of H] [391 HO] where fn,k = gn,]s., for all n and k. 
This results in = 0, which leaves only the new displacement interaction term in 
the Hamiltonian: Hj = Ho- In the variational case, however, the /„,k are left as free 
parameters, and minimization over an upper bound on the free energy, as described 
below, determines their values. The idea is that the optimization inherent to the 
variational approach allows us to minimize the effect of the interaction Hamiltonian 
Hi, given the transformation form. This is in done here in order to validate the use 
of perturbation series in various master equation approaches, which must in practice 
be truncated at some finite order. In general, the 'smaller' Hj, the more accurate the 
low-order dynamics are likely to be. 

Since, most of the time, there is no single parameter in Hi which determines 
exactly how small its effect is, one instead optimizes the variational transformation 
by minimizing the contribution of Hi to the free energy (the average energy of a 
thermal state of the system). As it is generally impossible to find an exact analytical 
expression for the free energy, it is the Feynman-Bogoliubov upper bound [19] that we 
shall minimize. The bound is given by 

^- ' ' + {H,}^„ + 0{{Hf}sV (8) 



trie 



where = tr (^Xe The true free energy A is related to this bound by the 

inequality A < Ab- Given that we want to end up with Hj small, it is reasonable to 
neglect the higher order terms in[8]as a first approximation. Furthermore, the interaction 
Hamiltonian in the transformed frame has been constructed such that the second term 
goes to zero, (Hj) = 0. Therefore, minimization amounts to maximizing the value 
of tr(e~^^"). Although, perhaps counter-intuitively, Hj now appears to be absent from 
Ab its influence is, in fact, still present implicitly in Hq. 

The transformed system Hamiltonian can be written as a function of the 
renormalization parameters {-R„,i3„}, therefore the minimization condition can be 
written: 

dAB_^dABdRr^^dABdB^^^ 

dfn,k dRn dfn,k 9/„,k 

which, after using the expressions for the renormalization parameters in [7] allows us to 
write fn,k = FniuJn,k)gn,k, with 



F„(c.„,k, {Rn, SJ) 2^^^^^ _ gj2^;^^h(/3a;„,k/2) ' 

In the continuum limit for the environment, the minimization procedure for an iV-site 
system amounts to solving the 2N coupled integral equations given by the definitions 
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of the renormalization parameters: 

Bn = exp -- du coth(/3w/2) , (11) 

\_ 2 Jo LO-^ 

Rn= doo -^^F^ioo, {Rn, i3„}) [F„(u;, i3J) - 2]. (12) 

Jo CO 

4. Master equation formulation in the variational frame 

As outlined above, a truncated perturbative expansion in the new interaction 
Hamiltonian Hj should, following the optimization procedure, be as accurate as possible, 
given the polaron form of the transformation and the minimization condition used. 
The next step, therefore, is to derive a master equation in the variational frame using 
standard techniques. By utilizing a projection operator V with the following action on 
the combined system-environment state: Vp = tr^(p) (g) pn, where p^ is an arbitrary 
reference state for the environment, we can separate out the reduced system dynamics. 
Here, we choose to derive a time-local or 'time-convolutionless' master equation due to 
the relative ease with which it can be solved numerically. In addition, we cut off the 
perturbation series at second order. The resulting master equation has the following 
general form in the interaction picture [28] : 

^ tlE {Vpm = tlE {}C2{t)Vp{t)} + tlE {X2(t)(l - P)p(O)} , (13) 

where /C2 and I2 are superoperators acting on Vp and {1—V)p respectively, which have 
been curtailed to second order in Hi. 

In the untransformed frame the separable initial state means that the second, 
inhomogeneous term in [13] disappears for the choice pr = p_b(0). In the variational 
frame this is no longer the case and the inhomogeneous term must be taken into 
account. However, for two-site systems the inhomogeneous term was seen to have 
only a small, transient effect on the dynamics at finite temperatures |44J for single 
site initial excitations. Therefore, we shall henceforth neglect it even in the transformed 
frame. This amounts to assuming that the environment relaxes into its displaced state 
instantaneously. One would expect this to be a good approximation when the typical 
environment timescales are much shorter than the transition timescales in the system. 

The remaining (homogeneous) term is written explicitly as 

tTE{IC2{t)Vp{t)} = - /ds tTE{[Hi{t),[Hi{s),rp{t)]]}- (14) 

By writing the interaction Hamiltonian in the form Hj = Y,i^=i Si ® Ei (with interac- 
tion picture counterpart Hi{t) = J2i^i Si{t) (g) Ei{t) ) we can rewrite the master equa- 
tion in terms of system operators Si and two-time environmental correlation functions 
Aij(t — s) = tiE {Ei(t)Ej{s)pR}. After moving back into the Schrodinger picture, the 
master equation takes the form: 
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= -i[Hs,ps{t)] 

Tds iA,,is){S,Sjis)ps{t)-Sjis)psit)S,} + h.c.). (15) 

The interaction Hamiltonian system operators can be split into three distinct groups in 
the following way: 

\n){n\ = l< i <N, 

\n){m\ + \m){n\ = S^^ N< i < |Ar(Ar + 1),(16) 

^ i\n){m\-i\m){n\ = lN{N + l)< i <N'^, 

which leads in turn to three varieties of non-zero time correlation function. The first 
type are due to the linear interaction term, Hi, and are therefore of the same form as 
those that appear in the standard weak coupling master equation: 

du Jn{uj) [1 - Fn{uj)f [cos{ujt) coth(/3a;/2) - isin(wt)], (17) 





where Fn{<jj) is the continuum version of the optimized function in [101 The second type 
come from the displacement operator interaction, if^i, and are the only type to appear 
in the fully displaced polaron master equation: 



+ exp [-S,^^iy{t)-5m,r;i!it)]-2 



^Z.pg{t) = \VnmV„BnBraBM CXp [^np^^^) + i^m, ' 5^pWI{t)] 



2 



exp [-dnpK'it) - (5^, - 5^,)<P^y{t)] , (18) 



where 



^iy{t) = / dcu -^^Fn{iof [cos(wt) coth(/3w/2) - isin(wt)], (19) 

JO OJ^ 

and the 6nm are Kronecker deltas. Finally, the third type appear in the more general 
variational master equation due to an overlap between the two types of interaction: 



with 



Ar„p(t) = 5r,^VnmBnBm<Pf{t), (20) 
'■^ , Jn(w) 



poo J 

{t) = / du ^^F„(a;) [1 - F„(a;)] [sin(a;t) coth(/3a;/2) + i cos(wt)] . 

JO U! 

(21) 

The dynamics calculated using [15] will be of the density matrix in the variationally 
transformed frame. In order to consider quantities in the lab frame, one must perform 
the inverse of the transformation in [H The site populations (diagonal elements of 
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the density matrix) are unchanged, since the operators commute with the 

transformation. However, in general, the inverse transformations of the coherences 
(off-diagonal elements of the density matrix) are much more difficult to calculate. In 
the case that the inhomogeneous term in [15] is ignored, one can make the approximation 
PU] pnmit) = BnBmPnmit), where p{s) is the system density matrix in the variational 
frame and p{s) is the system density matrix in the lab frame. This is equivalent to 
making a Born approximation in the transformed frame. 

5. Example system dynamics 
5.1. Three sites 

We would like to compare dynamics calculated in the variational frame, using [151 with 
other techniques. After two sites, the next simplest system with a Hamiltonian of the 
form of that in [1] has three sites and only nearest neighbour couplings. Figure [2] shows 
the dynamics for such a system in a variety of parameter regimes, calculated in the 
variational frame. For comparison, we have also plotted the dynamics calculated in 
the fully-displaced polaron frame as well as that calculated using the untransformed 
Hamiltonian (weak-coupling, or Redfield, approximation). The system is characterized 
by its on-site energies {En}, inter-site couplings V12 and V23, and spectral densities of 
the form 



Column (a) in figure [5] represents a regime where system frequencies (~ 20cm~^) are 
much smaller than the environment cutoff frequency (200cm~^). In this case, the 
fully-displaced polaron transformation is expected to do well [IS], and we would also 
expect the variational transformation to match it, as it indeed does. The untransformed 
dynamics fail to reach the correct steady state due to a reasonably large reorganization 
energy for the environment. 

Although the reorganization energy is the same in (b) as in (a), a weak coupling 
approximation (as performed in the master equation in the untransformed frame) in the 
former leads to unphysical results. This could be due to the fact that, by increasing 
one of the intersite couplings, one of the system eigenstate transitions has been put 
almost on resonance with the environment cutoff frequency. By better accounting for 
the effect of the environment, the variational and fuU-polaron transformations prevent 
the dynamics from suffering such a failure. 

The third and fourth columns, (c) and (d), show regimes for which neither the 
weak-coupling nor fuU-polaron transformation are ideally suited. System frequencies 
are comparable to environment frequencies, and the coupling to the environment is not 
small. The variational dynamics appear to interpolate between the two other results. It 
is clear that it agrees with the weak-coupling dynamics at short times, before settling 
on a different set of long time populations. In (c) these are much closer to the fuU- 
polaron results, but in (d) the comparatively large system energies prevent the polaron 




(22) 
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(a) (b) (c) (d) (e) 




5 2 4 0.5 0.2 0.4 0.2 0.4 



t/ps t/ps t/ps t/ps t/ps 

Figure 2. Population dynamics of a tliree-site system calculated using the variational 
(solid black), fuU-polaron (blue dotted) and weak-coupling (red dashed) master 
equations. The system parameters and spectral densities are as described in the 
text and in [211 with Wc = 200cm^^, Ei — E2 — 50cm~^, = Ocm~^ and 
T ~ 300 K in all cases. The remaining parameters for the individual subplots are: 
(a) Ai = A2 = eOcm^i, A3 = 120cm-\ V12 = V23 = 20cm-i; (b) Ai = A2 = 60cm-\ 
A3 = 120cm-i, V12 = lOOcm-i, V23 = 20cm-i; (c) Ai A2 60cm-i, A3 = 120cm-\ 
V12 = V"23 = lOOcm-i; (d) Ai = A2 = A3 = 180cm-i, V12 = 300cm-i, V23 = lOOcm-i; 
(e) Ai = A2 = A3 = 60cm- 1, V12 = 300cm-\ V23 = lOOcm^i. 

transformation from dealing correctly with the lower frequency parts of the environment. 
Additionally, in both cases, the dynamics in the fuU-polaron frame overestimate the 
damping of coherent oscillations. The variational transformation preserves coherence 
here precisely due to the fact that it optimizes the frequency dependence of the polaron 
transformation, as opposed to indiscriminately displacing every phonon mode by the 
full amount. 

Finally, in column (e) the reorganization energy is much lower than in column (d), 
but the system frequencies are still large. A weak coupling approximation is therefore 
valid in this case, and one would expect the dynamics in the untransformed frame to be 
more accurate than that in the (full) polaron frame. The red and black curves in this 
panel almost sit on top of each other, showing that the variational dynamics agree with 
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Figure 3. Variational (solid), weak-coupling (dashed) and fuU-polaron (dotted) 
dynamics for the populations of site 1 (black, red, blue) and site 3 (olive, orange, 
cyan) of the FMO complex. Although all seven sites were modelled, only the input (1) 
and output (3) site populations are shown for clarity. The system Hamiltonian used 
is the same as that in [50 and the spectral density is given in [231 Panels (a-c) have 
77=^, (d-f) have 77 1, (g-i) have rj = 2, and (j-1) have 77 = 4. The dynamics in the 
first column of plots were calculated at 5K, the second a 77 K, and the third at 300 K. 



the weak coupling results, and corroborates the fact that the variational transformation 
allows us to capture the dynamics across a broad range of coupling strengths. Note 
that, interestingly, both the weak-coupling and polaron approaches can overestimate 
the damping of coherence in comparison to the variational method, dependent on the 
parameter regime (cf. panels (a) and (e)). 

5.2. The Fenna- Matthews- Ols en complex 

We now analyse how the variational master equation performs for a larger system, 
namely the FMO complex. This is usually assumed to be a seven site system (although 
recent results suggest there is in fact an additional eighth site [8]) and thus has a much 
larger parameter space than the three site system in figure [2l Despite this, one can 
see from figure [3l which compares FMO dynamics across a range of reorganization 
energies and temperatures, that the variational transformation performs the same kind 
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of interpolation between weak-coupling and full-polaron dynamics as in the three-site 
case. The FMO system Hamiltonian used in this section is taken from [50] and the 
spectral density - the smooth part of that from - is of the form: 

J„(w) = 3.053 X 10"^ X r/^e-v^ + 1.908 x 10"^ x r/^e-v^, (23) 

with uji = 0.575cm~^ and U2 = 2cm^^. 

Figure |3] provides an clear example of the importance of the interplay between 
coherent and incoherent dynamics in excitonic energy transport. The biological purpose 
of the FMO complex is to transport excitations from site 1 (or sometimes site 6) to site 
3, from which the excitation is then removed [3]- It is therefore beneficial to have the 
population on site 3 build up as fast as possible. One can see from the figure (most 
clearly in the second column), that the optimum rate of transfer (panels (j), (h) and 
(i)) occurs in the variational theory when coupling to the environment is neither too 
strong nor too weak. That is, phonon-assisted transport is enhanced by the presence 
of some degree of coherence in the system. These optimal cases appear to lie in the 
intermediate region of parameter space, outside the remit of weak-coupling or polaron 
master equations, for which something like the variational approach is required. 

The full, seven-site dynamics for the cases where 77 = 1, at T = 77 K and T = 300 K, 
are shown in figure HI These plots correspond to the parameter regimes of the FMO 
dynamics in figure 2 of [50]. However, the exact calculations presented in that paper 
include a significant peak in the spectral density which, in any master equation approach, 
would ideally be treated separately from the rest of the environment in order to capture 
its effect on the dynamics non-perturbatively. Moreover, a variational polaron treatment 
of such a peak would likely cause there to be significant system-environment correlations 
in the transformed frame, which would not be taken into account without the inclusion of 
inhomogeneous terms in the master equation. That being said, the qualitative agreement 
of the variational dynamics shown in the figure with the results in [50] is surprisingly 
good. 

6. Discussion 

There are several advantages to using master equations over other approaches to open 
quantum systems dynamics. Primarily, these are the efficiency with which one can solve 
them, and the potential insights into underlying physics which they can give. We can 
see from the various terms in [15] exactly how the parameters in the Hamiltonian enter 
into the system dynamics and the relative magnitude of these terms can give us an 
idea of the parameter regime of a given system - in the sense of which quantities are 
contributing most to the dynamics. 

Moving into the variational frame prior to solving the master equation allows us 
to calculate sensible dynamics over a larger range of parameters compared to the more 
standard weak-coupling or full-polaron approaches. As can be seen from section [5l the 
variational master equation can capture the dynamics in both the weak and strong 



The multi-site variational polaron transformation 



13 



Pll P22 P33 PU P55 P66 P77 


\ (a) 


\ (b) 




:• • ^ 



I I. y- ' " I I I I 1 1, y- — » ■ - - I -t-. ........ 

0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 1 

t/ps t/ps 



Figure 4. Variational frame population dynamics for the FMO complex with the same 
parameters as figure 2 of [SOj, with the exception that only the smooth, non-peaked part 
of the spectral density was used. The dynamics in panel (a) were calculated at 300 K 
whilst those in (b) were calculated at 77K. Two different initial states: ps(0) = 
and ps(0) = 1 6) (6 1 were used for the upper and lower plots, respectively. 



coupling regimes, when the weak couphng and polaron master equations, respectively, 
are expected to do well. It is also able to bridge the gap between the two in intermediate 
regimes. 

Although the variational master equation is expected to work well over a wide 
range of parameters, and can in principle handle arbitrary spectral densities, it, like 
the full polaron transformation, works best in the scaling limit [45] - when important 
environmental frequencies (uc) are large compared to relevant system frequencies {Vnm)- 
The corresponding downside is that moving into the variational frame provides less of 
an advantage in terms of improving the accuracy of the dynamics when the typical 
environmental timescales are significantly longer than those of the system, and the 
coupling between the two is strong. The intuition for this is that the low frequency 
phonon modes are too 'sluggish' to keep up with the motion of the exciton as it moves 
through the system and do not, therefore, dress the system in the same way as higher 
frequency modes. They may still, however, have a profound impact on the system 
dynamics, which a transfomation of displacement form is unable to capture. 

One aspect of our method which might benefit from modification is the specific 



The multi-site variational polaron transformation 



14 



minimization condition used. Whilst we expect the first term of |S]to be a good metric 
for the size of the interaction, it does not directly correspond to the quantity which 
we expand perturbatively in the master equation. In fact, the next highest order term, 
{0{H])), is mathematically more similar [51], albeit much more complicated. As one 
smoothly varies the Hamiltonian parameters, the optimum transformation can jump, as 
different local minima become global minima. This effect has been studied for the case of 
two sites |17] , and it was found that the variational predictions are less accurate around 
the discontinuity. For multiple sites, the free energy landscape becomes more complex 
and more local minima emerge in parameter space, hence there is greater scope for 
this kind of jumping. Whilst not optimal, the transformations corresponding to such 
local minima are still likely to lead to better dynamics than those calculated in the 
untransformed frame. 

In summary, we have outlined a variational method for solving open quantum 
systems dynamics in molecular networks with local environments. The method is valid 
over a wide range of parameters and is efficient to compute. By moving into a reference 
frame in which system and environment are less strongly interacting, one is able to use 
a perturbative master equation to more accurately calculate dynamics. The approach 
surpasses both weak-coupling and fuU-polaron master equations in terms of breadth of 
applicability. It can be used to model interesting biological systems which sit in difficult 
intermediate coupling regimes, such as FMO, and allows for the systematic study of the 
effects of certain parameters on the dynamics. 

There is still ample room for improvement, and the general concept of redrawing 
the boundary between system and environment has far greater reach than the 
implementation presented in this paper. For instance, one could generalize the 
transformation to a larger class of Hamiltonians, say those whose environments couple 
to multiple sites, or one could augment the form of the transformation itself, perhaps 
by including squeezing in addition to displacement. The technique we have developed 
utilizes one of the most fundamental properties of quantum mechanics, namely the 
invariance of dynamical laws under unitary transformations, and gives insight into the 
important physical mechanisms underlying the evolution of open quantum systems. 
Whilst future master equation approaches may go beyond the polaron transformation, 
they are likely to benefit from a kind of variational minimization in the spirit of that 
which we have outlined here. 
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